//
// Created by Ivor_ on 2022/3/13.
//

#include "integral.h"

double igTrapezoidForPoints(int numberOfPoints, double** points) {
    double igv = 0.;
    for (int i = 0; i < numberOfPoints-1; ++i) {
        igv += fabs(points[i+1][0]-points[i][0]) * (points[i][1]+points[i+1][1]) / 2.;
    }
    return igv;
}

double igTrapezoidForFunction(int numberOfEqualParts, double start, double end, double* constantParameters, double (*pFunction)()) {
    double dx = fabs(end - start)/(1.*numberOfEqualParts), igv = 0.;
    if (end < start) {
        double temp = start;
        start = end;
        end = temp;
    }
    for (int i = 0; i < numberOfEqualParts; ++i) {
        igv += dx * (pFunction(start + 1.*i*dx, constantParameters) + pFunction(start + (1.+i)*dx, constantParameters)) / 2.;
    }
    return igv;
}

double igSimpsonForPoints(int numberOfPoints, double** points) {
    double igv = points[0][1] + points[numberOfPoints-1][1];
    for (int i = 1; i < numberOfPoints; ++i) {
        if (i%2) {
            igv += 4./3. * points[i][1] * fabs(points[i][0]-points[i-1][0]);
        } else {
            igv += 2./3. * points[i][1] * fabs(points[i][0]-points[i-1][0]);
        }
    }
    return igv;
}

double igSimpsonForFunction(int numberOfEqualParts, double start, double end, double* constantParameters, double (*pFunction)()) {
    double dx = fabs(end - start)/(1.* (numberOfEqualParts + (numberOfEqualParts%2)));
    double igv = pFunction(start, constantParameters) + pFunction(end, constantParameters);
    for (int i = 1; i < numberOfEqualParts; ++i) {
        if (i%2) {
            igv += 4.*pFunction(start > end ? end : start + 1.*i*dx, constantParameters);
        } else {
            igv += 2.*pFunction(start > end ? end : start + 1.*i*dx, constantParameters);
        }
    }
    igv *= dx/3.;
    return igv;
}

double igRombergForFunction(int order, double start, double end, double* constantParameters, double (*pFunction)()) {
    if (order <= 0) {
        order = 1;
    }
    if (end < start) {
        double temp = start;
        start = end;
        end = temp;
    }
    double* rombergArgs = (double*)calloc(order + 1, sizeof(double));
    rombergArgs[0] = .5 * (end - start) * (pFunction(start, constantParameters) + pFunction(end, constantParameters));
    int spaceSplit = 1;
    double hh = (end - start) / spaceSplit;
    for (int i = 1; i <= order; ++i) {
        double newPart = 0.;
        for (int j = 1; j <= spaceSplit; ++j) {
            newPart += pFunction(start + ((j << 1) - 1) * hh / 2., constantParameters);
        }
        rombergArgs[i] = rombergArgs[i-1] / 2. + hh * newPart / 2.;
        hh /= 2., spaceSplit <<= 1;
    }
    for (int i = 0; i < order; ++i) {
        for (int j = 0; j < order - i; ++j) {
            rombergArgs[j] = (pow(4., i + 1.) * rombergArgs[j+1] - rombergArgs[j]) / (pow(4., i + 1.) - 1.);
        }
    }
    double igv = rombergArgs[0];
    free(rombergArgs);
    return igv;
}

double igTrapezoidSelfAdaptation(double start, double end, double eps, double* constantParameters, double (*pFunction)()) {
    if (end < start) {
        double temp = start;
        start = end;
        end = temp;
    }
    double hh = end - start, igv[2] = {0., .5 * (end - start) * (pFunction(start, constantParameters) + pFunction(end, constantParameters))};
    int nn = 1;
    igv[0] = igv[1] + 99999.;
    while (fabs(igv[0] - igv[1]) > 3. * eps) {
        igv[0] = igv[1];
        double hn = 0.;
        for (int i = 0; i < nn; ++i) {
            hn += pFunction(start + (i + .5) * hh, constantParameters);
        }
        igv[1] = .5 * (igv[0] + hn * hh);
        nn <<= 1, hh /= 2.;
    }
    return igv[1];
}

double igRombergSelfAdaptation(double start, double end, double eps, double* constantParameters, double (*pFunction)()) {
    int order = 1, orderLimit = 1000;
    if (end < start) {
        double temp = start;
        start = end;
        end = temp;
    }
    double* rombergArgs = (double*)calloc(orderLimit + 1, sizeof(double));
    rombergArgs[0] = .5 * (end - start) * (pFunction(start, constantParameters) + pFunction(end, constantParameters));
    int spaceSplit = 1;
    double hh = (end - start) / spaceSplit;
    double newPart = 0.;
    for (int j = 1; j <= spaceSplit; ++j) {
        newPart += pFunction(start + ((j << 1) - 1) * hh / 2., constantParameters);
    }
    rombergArgs[1] = rombergArgs[0] / 2. + hh * newPart / 2.;
    hh /= 2., spaceSplit <<= 1;
    double igv[2] = {0., rombergArgs[0]};
    igv[0] = igv[1] + 99999.;
    while (fabs(igv[0] - igv[1]) > eps && order < orderLimit) {
        order++;
        newPart = 0.;
        for (int j = 1; j <= spaceSplit; ++j) {
            newPart += pFunction(start + ((j << 1) - 1) * hh / 2., constantParameters);
        }
        rombergArgs[order] = rombergArgs[order-1] / 2. + hh * newPart / 2.;
        hh /= 2., spaceSplit <<= 1;
        for (int j = order - 1; j >= 0; --j) {
            rombergArgs[j] = (pow(4., 1.*(order - j)) * rombergArgs[j + 1] - rombergArgs[j]) / (pow(4., 1.*(order - j)) - 1.);
        }
        igv[0] = igv[1], igv[1] = rombergArgs[0];
    }
    free(rombergArgs);
    return igv[1];
}

double igMcForFunction(int mcSteps, double start, double end, double* constantParameters, double (*pFunction)()) {
    double igv = 0.;
    for (int i = 0; i < mcSteps; ++i) {
        double loc = (start < end ? start : end) + fabs(end - start)*abs(rand())/RAND_MAX;
        igv += pFunction(loc, constantParameters)/mcSteps;
    }
    igv *= fabs(end - start);
    return igv;
}
